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We show how to systematically improve the Momentum Average (MA) approximation for the 
Green's function of a Holstein polaron, by systematically improving the accuracy of the self-energy 
diagrams in such a way that they can still all be summed efficiently. This allows us to fix some of 
the problems of the MA approximation, e.g. we now find the expected polaron+phonon continuum 
at the correct location, and a momentum-dependent self-energy. The quantitative agreement with 
numerical data is further improved, as expected since the number of exactly satisfied spectral weight 
sum rules is increased. The corrections are found to be larger in lower dimensional systems. 

PACS numbers: 71.38.-k, 72.10.Di, 63.20.Kr 
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I. INTRODUCTION 

One of the main challenges in condensed matter 
physics is to find simple yet accurate analytical approx- 
imations for systems in non-perturbative (strongly cou- 
pled) regimes. This is needed for interpretation of exper- 
iments and to help supplement the results available from 
numerical simulations of model Hamiltonians. 

In this general context, we have recently proposed the 
so-called Momentum Average (MA) approximation for 
the Green's function of a single Holstein polaron.^'— The 
essence of this approximation is to sum all diagrams con- 
tributing to the polaron self-energy, however each dia- 
gram is approximated to such a degree as to allow the 
analytical summation of the entire series. Specifically, 
each free propagator appearing in a self-energy diagram 
is replaced by its momentum average. The resulting 
MA self-energy is a trivial-to-evaluate continued fraction 
which gives remarkably accurate results over most of the 
parameter space, including intermediate electron-phonon 
coupling strengths where perturbational methods com- 
pletely fail to capture the correct physics. 

In Refs. [H and [3 we identified some of the reasons 
for this good agreement: first, the MA approximation 
becomes asymptotically exact for both very weak and 
very strong couplings. More importantly, the resulting 
MA spectral weight satisfies the first six spectral weight 
sum rules exactly, and remains highly accurate for all 
higher order sum rules. However, we also pointed out 
some shortcomings of the MA approximation: (i) it fails 
to predict the continuum that must appear at Eqs + f^, 
where Eqs is the polaron ground-state (GS) energy, and 
r2 is the frequency of the Einstein phonons (we set h = 1 
throughout this work). This continuum arises from states 
that have a phonon excited very far from where the po- 
laron is. As a result, their interactions are negligible and 
the energy of the system is simply the sum of the two. 
MA either predicts a wrong location for this continuum 
(at weak electron-phonon couplings) or no continuum at 
all in that range of energies (at intermediary and strong 
electron-phonon couplings). We noted in Ref. dthat nu- 



merical simulations show that there is very little spectral 
weight in this continuum, hence its absence or wrong po- 
sitioning does not significantly upset the agreement with 
the sum rules. Nevertheless, it would be reassuring to 
have an approximation that correctly predicts its exis- 
tence; (ii) the accuracy of the MA approximation wors- 
ens as ^ 0; (iii) the MA self-energy 'Sma{'^) is in- 
dependent of the momentum k of the electron. Given 
how featureless the Holstein model is (electron-phonon 
coupling and phonon frequency are both constants) one 
may expect a rather weak momentum dependence of the 
self-energy, however it is certainly not entirely absent. 

In this article, we show how to systematically improve 
the MA approximation, generating a hierarchy of approx- 
imations that we call MA^") (the original MA is MA^") 
in this notation). As explained below, the idea is to sys- 
tematically improve the accuracy of the "simplified" self- 
energy diagrams. The results become more and more 
accurate as n increases - for example, while the MA spec- 
tral weight satisfies only the first 6 sum rules exactly, this 
improves to 8 and 10 exact sum rules respectively for 
MA(i) and MA^^) spectral weights. While the numeri- 
cal effort also increases, it is still trivial for the n = 1 
and n — 2 levels that we discuss explicitly here. Level 
MA^^) already solves the continuum problem, while all 
levels with n > 2 produce momentum-dependent self- 
energies. The accuracy in the limit f2 — > is also shown 
to improve significantly with increasing n. In effect, for 
a slightly increased numerical effort, MA*^^^ solves all the 
known problems of the MA*^°^ approximation. 

The work is organized as follows: in Sec. II we briefiy 
review the MA(o) approximation, presenting a new ar- 
gument to explain its accuracy. In Sec. Ill we describe 
the systematic approach to obtain the improved versions 
MA("), n > 1, and give explicit formulae for the self- 
energies corresponding to the MA^^^ and MA*^^) approxi- 
mations. In Sec. IV we compare the predictions of these 
approximations against numerical simulations, to gauge 
the improved accuracy as n increases. Spectral sum rules, 
as well as variational arguments, will also be used to ex- 
plain the systematic improvement of accuracy with in- 
creasing n. Finally, Sec. V contains our conclusions. 
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II. BRIEF REVIEW OF THE MA^"^ 
APPROXIMATION 

The Holstein model is the simplest lattice model that 
includes electron-phonon coupling. Its Hamiltonian is:^ 



kC^Ck 



k,q 



b-r 



The first term is the kinetic energy of the electron, with 
cjj. and Ck being the electron creation and annihilation 
operators. For the single electron (polaron) problem of 
interest to us, the spin of the electron is irrelevant and wc 
suppress its index. Ck is the free-particle dispersion. In 
all results shown here, we assume nearest-neighbor hop- 
ping on a d-dimensional simple cubic lattice of constant 
a (we set a = 1) with a total of sites, and with periodic 
boundary conditions. In this case 



Ek = — 2ty~^ cos(fcia) 



but our results are valid for any other dispersion. The 
second term describes a branch of optical phonons of 



energy ft. &^ and 6q are the phonon creation and an- 
nihilation operators. The last term is the on-site lin- 
ear electron-phonon coupling V|:i_ph — gJ2i + ^«)' 
written in k-space. All sums over momenta are over the 
first Brillouin zone, namely — tt < ki < tt, i = l,d. 

The quantity of interest to us is the (retarded) single 
polaron Green's function, defined asi^ 



GiKu) = (0|ckG(u;)4|0) 



1 



Ck- 



n + ir] 



4lO) (1) 



where |0) is the vacuum Ck|0) = &q|0) = 0, and 77 > is 
infinitesimally small. 

As described in detail in Ref. 0, using repeatedly 
Dyson's identity G{uj) — Go(cj) -I- G{uj)Vc\-phGQ{uj), 
where Vd-ph = Ti. — TLq is the electron-phonon interac- 
tion potential and Go(a') = [lo — Ho + «'7]~^, we generate 
the infinite hierarchy of equations of motion whose exact 
solution is the desired Green's function: 



G(k,c^) = Go(k,c^) 



and for n > 1, 



qi 



, (2) 



J 



i^„(k,qi, . . . ,q„,tj) = -j=Go(}i~ qr, uj-nn) 



.1=1 



qn+i 



(3) 



Here, — 1* is the total momentum carried by 

phonons, Go(k, w) = (w — Ck + *?y)^^ is the free elec- 
tron Green's function, and wc introduced the generalized 
Green's functions 



FnO^, qi, 



\c^.G{u;)ci^^X^■■■^i^^)■ 



These equations can be recast in a more convenient 
form after observing that if we treat Eqs. ([3]) as an 
inhomogeneous system of linear equations in unknowns 
Fi,F2,..., then the only inhomogeneous term appears 
in the first equation and is proportional to Fo(k, w) = 
G(k, w). It follows that all generalized Green's functions 
Fi,F2,... must be proportional to G(k, w). As a result 
we introduce the more convenient variables: 



/«(k, qi, 



7V35"i^„(k, qi, . . . ,q„,w) 



G(k,w) 

In terms of these, Eq. ([2]) becomes: 



(4) 



G(k,t^) = Go(k,tj) 



N ^ 
qi 



/i(k,qi,c^)G(k,c^) 



SO that the exact self-energy is: 



E(k,c.) = l^/i(k,qi,o;), 
qi 



(5) 



giving the standard solution: 



G(k,w) 



w — ek — S(k, a;) + ir\ 



(6) 



To find /i (k, qi , , we must solve the infinite sys- 
tem of coupled equations that result from Eqs. ([3]). 
For later convenience, we write the first few equa- 
tions explicitly here, using the short-hand notation 
/„(k, qi, . . . ,q„,Lj) = /„(qi, . . . ,q„) (i.e., the depen- 
dence of k and lo of these functions is implicitly assumed 
from now on). Then: 
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and for all n > 3, 



/i(qi) = Go(k - qi,w - Q) 



^5I-^2(qi,q2) 
q2 



/2(qi,q2) = Go(k - qi - q2,t^ - 2n) 



f [/i(qi) + /i(q2)] + ^5I-^3(qi,q2,q3) 



/„(qi, . . . ,q„) = Go(k-^q,,w-n^2) 



4=1 



(7) 
(8) 



" 1 



qn+1 



(9) 



Clearly, all the dependence on free propagators of the 
self-energy comes from the free propagator prefactors on 
the right-hand side of these equations. 

As already stated, MA^^^ consists in replacing all free 
propagators in all self-energy diagrams by their momen- 
tum average: 



(10) 



Obviously, this corresponds to replacing the free propa- 
gator pre- factors on the right-hand side of Eqs. Q-© by 
the corresponding go(<^ ~ nil). In this case, it is straight- 
forward to see that the resulting solutions, denoted 
are functions of lu only (all dependence on phonon mo- 
menta disappears at this level of approximation). The re- 
sulting equations /;[°''(a;) = go (^ — ^^) g"^ + f2^\^) and 
for n > 2, fL°\Lu) = goito-nn) [ng^fi'Uu;) + fi%{u;) 

are solved in terms of continued fractions^ (also see Ap 
pendixEl to findi^ 



where we define the infinite continued fractions: 

ngo(LU — nQ) 



1 - .9^.90 (t^ - nn)A,i+i{uj) 
ngo{LU — nQ) 



(11) 



(12) 



{n + l)g go{uj — nQ)go(LU — [n + l)fi) 



1 - . . . 

Results based on this MA'^*'^ approximation have been 
analyzed in detail in Ref. 2. 

Before discussing how to systematically improve this 
approximation, it is worth pointing out an alternative ex- 
planation for its good accuracy.'^ This comes from realiz- 
ing that in real space, the meaning of the MA^^^ approxi- 
mation is that it replaces all free-propagators Go{i,j,Lu — 
jifl) appearing in all self-energy diagrams, by dijga{uj — 
nfl). For example, consider the real-space, second-order 
Green's function diagram depicted in Fig. [T] It has the 




FIG. 1: A second order diagram contribution to G{i,j,uj) 



exact value X^iiji G'o(i,«i,w)S2,c(«i, ji,t^)Go(ji, j,tj), 
where the contribution to self-energy from the second- 
order crossed self-energy diagram is: I]2.c(*i, ji, ^^) = 
9^Go{ii,ji,uj - il)Go{ji,ii,uj - 2n)Go{ii,ji,uj - n). 
Within MA^''^ I]2,c(*i, ji, w) is approximated as 



X. . y(0) 



ico) 



IS 



il)go{Lu 



2fl)go{u! — n). Inserting this into the Green's function di- 
agram removes one of the sums, and after Fourier trans- 
forming we find that the contribution of this diagram to 
G(k,w) is Go(k,w)S^°^(w)Go(k,w). The same holds for 
all higher order diagrams. Summing all of them, we see 
that the free propagators in the proper self-energy parts 
have indeed been replaced by their momentum averages. 

The reason why it is a good zero-order approximation 
to keep only the diagonal (in real space) contributions 
Go(i,j,i^ — riQ) Sijgo{uj — nft) is straightforward to 
understand, at least for low energies ui ~ Egs- Because 
of interactions, Eos < —2dt (the polaron ground-state 
is below the free electron continuum). It follows that for 
uj ^ Egs, the free propagators Go(«, j, uj—nil) are needed 
at energies below the free electron continuum, and the 
larger n is the further below the band-edge these energies 
are. However, it is well-known that for uj < —2dt, the 
free propagator decreases exponentially with increasing 
distance \i — j\. The reason is that there are no free- 
electron eigenstates outside the free-electron band, and 
the electron has to tunnel from one site to another. For 
example, in one dimension^ 

Go(z,j,c.) = e*'l^-^"lGo(*,z,c^) = e^^-«l'-^lgo(t^) 

where fco is the first quadrant solution of the equation 
Lo = —2tcosko. For w < —2t, ko = iK, where k > 
increases as cu decreases. We also used the fact that in 
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terms of its Fourier transform: 



It follows that at low energies, MA^"^ keeps the largest 
contributions to the self-energy (from real-space diagonal 
terms), while ignoring exponentially small contributions 
coming from off-diagonal terms. This is expected to be- 
come more and more accurate for higher order diagrams 
with many phonons: the larger n is, the faster the expo- 
nential decay with distance of Go{i,j,LO — nil) becomes. 

Based on these arguments one expects at least low- 
energy properties to be well described by MA, at least if 
51 is not too small. Together with the good agreement-!^'^ 
with the spectral weight sum rules, this leads to the con- 
clusion that the self-energy and the Green's function at 
all energies should indeed be quite accurate within MA. 
Comparison with numerics validates thisi^ 

III. HIGHER LEVELS MA("\n > 1 

The above arguments also suggest a systematic way to- 
wards improving the MA^"-* approximation. The biggest 
error at low energies is due to the momentum average 
of the propagators of energy uj — fl. This is the en- 
ergy closest to the free electron continuum and therefore 
these propagators have the slowest decay in real space. 
The next slowest decay is for the propagators of energy 
u! — 251, etc. If one could selectively keep these propa- 
gators exactly while momentum-averaging the ones with 
more phonons (lower energy, faster decay) this should im- 
prove the accuracy of the approximation at low energies. 
In fact, all sum rules would also be improved (individual 



diagrams are more accurate) therefore one would expect 
an improvement at all energies. This is precisely what 
the higher levels MA("\ n> 1, achieve. 

We define MA^"^ as being the approximation where in 
all self-energy diagrams, all free propagators with energy 
uj — mil, where m < n, are kept exactly, while the ones 
with more than n phonons are momentum averaged. 

In terms of the equations-of-motion that need to be 
solved, Eqs. ©-(O, achieving this is straightforward, 
since a propagator of energy uj — niVL appears only once, 
in the right-hand side pre-factor of the equation for 
/,„(qi, . . . , q„j). It follows that if we keep the first n 
of Eqs. ([7])- ([5]) as they are, and approximate the equa- 
tions for fn+i, fn+2, • ■ ■ by momentum-averaging the free 
propagator appearing in the right-hand side pre-factor 
Go(k — qr,!^ — mfl) — > goii^ — mfl) if m > n -|- 1, we 
achieve our goal - provided that we can find the solution 
of the resulting infinite system of coupled equations. 

We now derive explicitly the solutions for n = 1 and 
n — 2 levels. 

A. MA(^) level 



In this case, the equations to be solved are: 
EM^(i,(k,c.) = ^5]/l^)(qi) 



(13) 



qi 



where 



/f^(qi) = Go(k-qi,c^-51) 
and for all n > 2, 



(14) 



N 



/P(qi, ■ • ■ ,q«) = ffo(w - nfi) 



i=l 



qn+i 



(15) 



As before, the dependence on k, lu is implicitly assumed 
everywhere. We use the upper labels (1) because these 
are the approximative solutions corresponding to MA'^^^ 
The solution of this infinite set of recurrent equations 
is discussed in Appendix IB] The end result is: 

1 - 9 {MK^) - Ax(uj - S1)J 

where [see Eq. (fTTj) ]: 

UJ — UJ — ^ — A\(uj — 51) = cij — 51 — Ema('^ — f^)- (17) 

The continued fractions Ax(lo — 51), A^ito) are defined in 
Eq. p2p . This expression is slightly more complicated 



than Sj(,f^(o) (tij), since it involves two different continued 
fractions, however it is still very trivial to compute. 

Note that based on this and other results de- 
rived in Appendix [B] we can now calculate the 
expressions for the generalized Green's functions 
/l^^(k, qi, . . . ,q„,w) - F„(k, qi, . . . ,q„,w). These will 
be more accurate than the values obtained within the 
MA*^") approximation, where none of the expres- 
sion had any momentum dependence. These generalized 
Green's functions contain further information about the 
polaron, for example regarding the phonon statistics. 

As can be seen from Eq. ((T^ . the self-energy is still 
momentum independent at this level. However, it is 
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clear that this is because the Holstein model is so fea- 
tureless. If, e.g. the coupling was dependent on the 
phonon momentum, then the first self-energy diagram 
|g(q)pGo(k — q, would be k dependent, and 
so would "Ej^^j^ii) (k, Lu) (this diagram is exact at MA'^"'^^ 
level). Indeed, work in progress on generalizing this ap- 
proach to models with .g(q) coupling verifies this. 

It is also clear that even for the Holstein model, all ex- 
pressions (k, Lu) with n > 2 will have momentum 
dependence, since Holstein self-energy diagrams of sec- 
ond order are momentum dependent. We demonstrate 



J 



below that this is indeed the case. 

Finally, the MA approximation becomes exact in the 
limit g — > and t — > 0. The first limit is trivial, since 
S ^ 0. The second is due to the fact that ii t = Q 
then free propagators Go(k, w) are in fact independent 
of k, and thus the momentum averages become irrele- 
vant. Clearly, the same must hold true for all higher 
level MA^") approximations. Indeed, one can verify di- 
rectly that ifgo{oj) = {Lj + ir])~^ (corresponding tot — 0), 
then the expressions in Eqs. and pT|) are equal. The 
same is true for the MA^^^ results we present below. 



B. MA(2) level 



In this case, the equations to be solved are [compare to the exact Eqs. ([T])-®]: 

EMA(^)(k,c.) = l^/p)(qi) 



where 



/f^(qi) = Go(k-qi,c.-r!) 



/f ^(qi, qa) = Go(k - qi - qa, - 217) 



5' + ^E/i'^(qi'q^) 

/f^(qi) + /f^(q2)l + ^E/f (qi'^^^qa) 



(18) 

(19) 
(20) 



and for all n > 3, 



/l^^(qi, . . . ,q«) = goiuJ - nfl) 



n ^ 

9^ XI ^"-1 ■ ■ ' ' • • ■) + ^ XI (qi ' • • ■ ' qn-t-i) 



i=l 



(21) 



Dependence on k, oj is again implicitly assumed. 

This can be reduced to a closed system of equations in 
terms of only f{ (qi) and /2 (qi,q2), after solving for 

:^Eq3 /3^^(qi'q2,q3) from Eqs. Hi]). The details are 
provided in Appendix [Cl 

We use the short-hand notation: 

Ai = Ai{uj- 2Vl)-A2 = A2{uj- n)- As = A3(w), 

where the continued fractions are defined in Eq. 

We also define various momentum averages (dependence 

on k, is implicit): 



•^1 - ^MA(2) - 1^ /i^^(qi)i 

qi 

1 ^ J-?^. , 2g^go{Cj)Ti 



'5/2(qi) = ^E/2'^(qi'q2)-^2. 

q2 

The link between J-i and J-'2 is proved in Eq. (|C4 



In terms of these, the closed system of equations to be 
solved becomes (see Appendix [Cl for more details): 

/f ^(qi) = Go(k - qi, - f7) [g' + 5/2 (qi) + , 



5/2 (qi) = .9^.90 (^) /f ^(qi) + (^2 - A,)S h{qi) - 2Tr 



' N 



^Go(k-qi-q2,(:>) [/(')(q2) + (A2-Ai)(5/2(q2) 
q2 

These equations can be solved in a variety of ways. 
We present here the most efficient solution that we have 
found, and then comment briefly on other possible so- 
lutions. First, given the form of these equations, it is 
advantageous to introduce the new unknown 



fi'\q) + iA2-Ai)6f2{q). 



(22) 



Consider its Fourier transform at various lattice sites Ri, 
namely a:(i) = jfJ2q^^'^^'^i- First, observe that a;(0) = 
ir Sq 2^q = -^^i = ^MA(^) I by definition. 
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As shown in Appendix [Cl the set of two closed equa- 
tions above can be rewritten as an inhomogeneous equa- 
tion involving x{i) (at all lattice sites): 

Y^M,j{k,uj)x{i) = e'^'-^^g^GoH,^) (23) 



where 



1 - g'^go{u){A2 - Ai) 



Go{i,uj) — e''^'^*Go(k, cj), and the expression of 

the matrix elements My(k,a;) is given in Eqs. (|C5IIC7[) . 
They have simple expressions, involving only (the same) 
three continued fractions Ai, A2, A3 as well as various 
Go{i,oj) values, therefore they can be calculated easily. 

Because at low energies the free propagators in real 
space decay exponentially, one expects that x{i) de- 
creases fast with increasing distance R^. Alternatively, 
consider, for instance, the /j (q) contribution to Xq. 
When Fourier transformed, the initial state cj^_qfcq|0) 

goes into Cjbj_^^\0), i.e. the phonon is further and fur- 
ther apart from the electron. Similar interpretation can 
be given for the second contribution to x{i). One expects 
the amplitudes for such processes to decay with i. 

As a result, we can truncate the system of equations 
assuming that x{i) = for all larger than a cut- 
off. We vary this cutoff to insure that convergence has 
indeed been achieved. In this formulation, we find that 
convergence is reached extremely fast, typically for a cut- 
off distance of order 5a (see results section). In other 
words, to obtain E^^^(2) (k, ijj) — x(0), in ID we typi- 
cally have to solve a system of 11 or so inhomogeneous 
equations, which can be done very efficiently. Higher di- 
mensions imply larger systems, but overall the numerical 
task is still trivial and results can be obtained very fast 
and with little computational resources. 

It is important to emphasize that such low cutoffs are 

not inherent in the problem. In fact, one could also solve 

(2) 

these equations, for instance, by ehminating /| "^(q) to 
get an equation only in terms of i5/2(qi) and ^ S. If 
one Fourier transforms this, it turns out that cutoffs as 
large as 100a are needed before convergence is achieved, 
and this is especially so for the polaron+one-phonon 
continuum (the bound polaron states converge quickly). 
This is not surprising, since states in the polaron-|-one- 
phonon continuum do have a free phonon, i.e. one that 
could be infinitely far from the polaron. In reality, we ex- 
pect that if we have a finite but large enough system, all 
quantities will eventually converge to their bulk values. 
In particular, here this suggests that one needs to allow 
the free phonon to be hundreds of sites away from the po- 
laron before convergence for the continuum is achieved. 

The much faster convergence for the formulation of Eq. 
(1^ is due in part to the particular choice of variable 
Xq. Even more important is the infinite summation of 
diagrams performed when the new frequency uj appears 



(see Appendix [C]). Without this, the convergence for 
continuum energies remains slow and large cutoffs are 
needed. Examples are discussed in the results section. 

Of course, one could also attempt to solve these equa- 
tions directly in the k-space. Without having tried it, 
we believe this to be an inefficient approach. The goal is 
to find J^i = S, i.e. an average over the Brillouin zone. 
Within MA(i), f[^\q) ~ G'o(k - q, w) (it is a constant 
in MA). Presuming that /f (q) is not too different, it 
is clear that these functions are of comparable size ev- 
erywhere in the Brillouin zone, and therefore one should 
sample many points in the Brillouin zone to obtain an 
accurate average. 

Finally, going back to Eq. , we would like to point 
out that if we set the cutoff at zero distance, i.e. use 
Mqox{0) — g'^g{){uj), we obtain an analytical, momentum- 
independent approximation to the true S^^(2) (k, w): 



where 



3^90 (w) 



2go(ci)5o('I') [-^^^ - 



, (24) 



aijiuj) = 1 - g^go{Q){Ai - Aj). 

One can think of this as the variant of MA*^^^ that keeps 
some of the free propagators of energy w — 2i7 exactly 
(typically those appearing in non-crossed diagrams) , but 
averages over those that give momentum dependence to 
the self-energy. The self-energy 'E]^j^^(2) (w) is more accu- 
rate than {to) but less accurate than E^^(2) (k, uj). 
Such "zero-cutoff" analytical approximations can be ob- 
tained for higher levels of MA*^"^ quite easily. The full 
MA(") for n > 2 can also be done. The reduction to a 
closed system of n coupled equations is always straight- 
forward. Its solution, however, becomes more computa- 
tionally involved as n increases, and leads to gradually 
smaller improvements in the accuracy. 



IV. RESULTS 

A detailed comparison of the predictions of the MA ap- 
proximation vs. numerical simulations is available in Ref. 
(2|. Instead of another comprehensive investigation, here 
we will focus on several properties where the higher level 
MA*^") approximations show a significant improvement 
over MA results. The way of extracting quantities of in- 
terest from the Green's function, e.g. ground-state (GS) 
energies, quasiparticle (qp) weights, effective masses, av- 
erage number of phonons in the polaron cloud, etc., are 
described in detail in Ref. 2. Throughout we use 



A = 



2dtn 



(25) 



as the effective coupling strength, d being the dimension 
of the lattice. All energies are measured in units of t. 
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FIG. 2: (color online) (a) Ground state energy and (b) 
Ground state quasiparticle weight, as a function of the ef- 
fective coupling A, for d = l,t = l,n/t = 0.1. The QMC 
results are from Ref. [l^. 



A. Ground-state properties 

The ground-state energies predicted by MA are quite 
accurate for a large range of parameters. The accuracy is 
known to worsen as fi/f — > 0, however even for il/t = 0.1 
the MA energies have less than 5% relative error. On the 
other hand, the GS qp weight for this low fl/t ratio is 
quite wrong for intermediary couplings A ~ 1, although 
it does become asymptotically exact, as expected. The 
comparison with Quantum Monte Carlo (QMC) results 
is shown in Fig. [2l where we also show the MA^^^ and 
MA^^) predictions. 

The accuracy is improved significantly for the higher 
MA levels. Improvements are observed for all other sets 
of parameters (not shown) analyzed in Ref. d, however 
for higher Q/t ratios MA is much more accurate to begin 
with, so the supplementary improvements due to MA'^' 
and MA^^^ are comparatively smaller. As discussed, this 
systematic improvement is expected since all self-energy 
diagrams become more and more accurate. This is re- 
flected in the sum rules for spectral weight, which are 
also systematically improved (the link between diagrams 
and sum rules is discussed at length in Ref. 0). While 
MA satisfies the first 6 sum rules exactly, MA^^^ satis- 
fies the first 8 sum rules exactly, MA^^^ satisfies the first 
10 sum rules exactly, etc. Of course, the accuracy of all 
higher order sum rules is also systematically improved by 
the use of more accurate expressions for the diagrams. 

This systematic improvement can also be understood 
in variational terms. As pointed out in Ref. ^ MA is 
equivalent to a variational approach where the eigenfunc- 
tions are built within a restricted Hilbert space that al- 
lows phonons only at the site where the electron is. In 
other words, the real-space counterpart of the general- 
ized Green's functions F„, namely Fn{i; j, ji, . . . , j„; = 



{0\cM^)c]bl---bl\0) -> F^°\t;j;u;)Ul=iS,,,, at the 
MA(°) level. This is equivalent to asking that the single- 
polaron eigenstates have non-zero overlap only with basis 
states of the general type c](&J)"|0), Vj, n. It is straight- 
forward to verify that with this restriction, the resulting 
equations for F„ lead to the MA'^'^^ self-energy. 

This observation immediately explains the absence of 
the polaron-|-one-phonon continuum, since within this re- 
stricted Hilbert space it is not allowed to have a phonon 
far from where the electron is. To compensate for the 
missing continuum's weight and satisfy the sum rules, 
the GS qp weight is increased within MA, as seen in Fig. 
[2Jb) and many other examples discussed in Ref. 

In this variational interpretation, MA^^^ almost corre- 
sponds to using a restricted Hilbert space enlarged by 
basis states of the form cj(6])"'&j, |0), Vj, i.e. it 

also includes states where one phonon could exist any- 
where with respect to the electron. The equivalence is 
not exact, because the resulting equation for F2 is not 
the same in the two cases (all other equations for all 
other Fn with n 7^ 2 are the same). More precisely, one 
should think of MA^^^^ as a variational method also ac- 
companied by a change in the Hamiltonian if and only if 
acting on electron-|-two-phonon states. Similarly, MA*^^^ 
corresponds to using a restricted Hilbert space spanned 
by basis states that allow any number of phonons on the 
electron site plus up to two phonons anywhere else in the 
system, accompanied by a change of the Hamiltonian if 
and only if acting on electron-|-three-phonon states (now 
the equation for F3 is not quite the same as in MA*^^^), 
etc. Note that one could also define improvements to 
MA based on the variational equations for F2 (instead of 
that resulting from MA^^^), F3 (instead of MA^^)), etc. 
However, these equations are more involved than the cor- 
responding MA^") type equations, making the solution of 
the system of coupled equations more difficult. 

This systematic increase in the size of the variational 
Hilbert space is another way to explain the gradual im- 
provement of the GS energy. Also, it is now clear that 
the polaron + one-phonon continuum should appear at 
level MA*^^) (see below). As a result, the GS qp weight 
no longer needs to account for it and it decreases, im- 
proving the agreement with QMG results as shown in 
Fig. Hfb). We will return to this issue when we investi- 
gate spectral weights. For the time being, we note that 
in the adiabatic limit —^ 0, many phonons can be 
created in the system at low energetic cost. In the in- 
termediary region A ~ 1 where the polaron cloud is still 
relatively large, one expects many of these phonons to 
be relatively far from the polaron and therefore a high 
order n would be required in order to accurately describe 
them within this approach. As a result, it is expected 
that the strongly adiabatic regime will not be quanti- 
tatively well described for A ~ 1 by the low-level MA 
approximations, even though the qualitative behavior is 
correctly captured. Of course, this limit can be investi- 
gated by other means, such as in Ref. [Tj and references 
therein. However, for most of the parameter space, i.e. 
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FIG. 3: (color online) (a) and (b) Polaron dispersion S^; (c) 
and (d) qp weight Z^, and (e) and (f) average number of 
phonons Nph{k) vs. fc, in d = 1, for Q = 0.5t and A — 0.25 
in (a),(c)(e), respectively A = 1.00 in (b),(d),(f). The QMC 
results are from Ref. Hsl. 



any Q/{dt) > 0.1 or so and any coupling A, the MA set of 
approximations give very easy to evaluate yet remarkably 
accurate results. 



B. Polaron band 

We can also track how the lowest eigenstate of momen- 
tum k, and its various properties, change with various 
parameters. Results are shown in Fig. [3] for the energy 
i?k, IP weight Zk and average number of phonons A^ph(k) 
for ID and two couplings, A = 0.25 and A = 1.00. We 
found similar improvements in 2D cases. Clearly, MA*^^' 
leads to an obvious improvement, even though for this 
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FIG. 4: (color online) (a) Spectral weight A{k = 0,u}) vs. uj 
in ID for f = 1,0 = 0.5, A = 0.6, r? = 0.01, in MA, MA'^' 
and MA''^-' (curves are shifted for clarity); (b) Polaron+one- 
phonon continuum convergence with cutoff within MA'^' . Re- 
sults for cutoffs of 0, 1, 3, 5 and 10 are shown (the last two 
are almost identical). For comparison, the MA'^' continuum 
is also shown (black full line); (c) same as in (b), but for an 
inefficient computation scheme. 



value of fl/t = 0.5, MA itself is quite accurate already. 
It should be noted that A ~ 1 is where the MA accuracy 
is generally expected to be at its worst. In particular, for 
the weak-coupling value A — 0.25, we see that MA overes- 
timates the distance to the continuum, i.e. i?^ — Eq > Q. 
This should be fl, but it is larger for MA because the 
polaron-l-one-phonon continuum is not predicted at the 
correct energy (For A = 1, there is a second bound state 
between the one shown and the continuum, therefore the 
bandwidth is much less than fi). For MA^^^ and MA^^^ 
this problem is indeed fixed, and the polaron dispersion 
width (at weak couplings) is Q. All other quantities are 
also clearly more accurate. 



C. Polaron+one-phonon continuum and higher 
energy states 

In order to understand the effects on higher-energy 
states, we study the spectral weight A{]s.,uj) = 
— ilmG'(k, Lu). As is well known, this is finite only at en- 
ergies Lo where eigenstates of momentum k exist. For dis- 
crete (bound) states the spectral weight is a Lorentzian 
of width 77 and height proportional to the qp weight. In a 
continuum, the lifetime is determined by Iml](k, cj) and 
is independent of 77 if 77 is chosen small enough. 

In Fig. m^a) we show results for the ID spectral weight 
A{k = 0, Lu) vs. LO for a relatively weak coupling A = 0.6. 
The MA spectral weight shows two discrete states at low 
energies, and a continuum starting for lo > —l.bt. Within 
MA'^^^ the second peak spreads into a continuum whose 
lower edge is at roughly above the energy of the GS 
peak. In fact, since go(w) acquires an imaginary part 
when —2dt < lo < 2dt, from Eq. (fTB|) it follows that the 
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FIG. 5: (color online) (a) Spectral weight ^1(0, ui) vs. uj in ID 
for t = 1, n = 0.5, A = 1.2, ri = 0.01, in MA, MA^^' and MA'^' 
(curves shifted for clarity); (b) lnj4(0,tj) vs. lj for MA^^' and 
r] = 10~^, 10~^, lO"**. Other parameters are as in (a). 




FIG. 6: (color online) (a) Spectral weight A{0,uj) vs. uj in ID 
for t = l,Q = 0.5, A = 1.8, ri = 0.01, in MA, MA^^' and MA'^' 
(curves shifted for clarity); (b) lnj4(0,tj) vs. lj for MA(2' and 
rj = 10~^, 10"'^, 10~^. Other parameters are as in (a). 



MA*^^^ continuum appears when uj > — 2t ^ a; > E^g + 
ri, where E^g is the MA prediction for the GS energy. 
Since the MA^^) GS energy E^^l < e''°1, it foUows that 
the gap is in fact slightly larger than fi. In MA^^^ the 
weight is redistributed within this continuum, and its 
lower band-edge is at f2 above the ground-state energy, 
within numerical precision. 

The convergence of MA^^^ with the cutoff value used 
to truncate Eqs. (|23p is shown in Fig. [UJb). For a cut- 
off value of we obtain the momentum-independent self- 
energy of Eq. (|24p , which gives a continuum with a shape 
rather similar to that predicted by MA^^^ As the cutoff 
value is increased, weight shifts towards the lower band- 
edge. Convergence is reached very quickly, with little dif- 
ference visible between results corresponding to a cutoff 
of 5 or 10 (these imply solving an inhomogeneous system 
of 11, respectively 21 equations in Eqs. (|23|) ). Other 
solutions of the coupled equations (briefly discussed in 
the previous section) converge much more slowly. In Fig. 
mjc) we show results for 3 cutoffs for such an alternative 
scheme. Even for a cutoff as large as 100, one can still 
see small oscillations, very reminiscent of finite-size ef- 
fects. The finding of an efficient solution for the MA^^' 
self-energy is thus quite important. 

Similar behavior is observed for higher couplings, as 
seen in Figs. [5fa) and[Bfa) for intermediate, respectively 
strong couplings A = 1.2 and 1.8. In both cases there 
are now two bound, discrete states below the continuum 
starting at Eqs -I- as expected.™ The weight of the 
polaron-|-one-phonon continuum decreases very fast, so 
that for A = 1.8 it is barely visible just above the sec- 
ond peak. Fig. EUb) shows it clearly, on a logarithmic 
scale. Interestingly, it is not only the height of this fea- 
ture that is much smaller as A increases, but its width as 
well. In fact, scaling vs. rj in Fig. [G^b) shows that this is 
more like a Lorentzian, i.e. a single bound state, and not 
a finite-width continuum as was the case for lower cou- 
plings (see Fig. (Hb)). This is in fact reasonable, since 



at such large couplings the lowest energy polaron state is 
basically dispersionless (the effective mass is already con- 
siderable and the polaron is well into the small-polaron 
regime). Since the width of the continuum is due to the 
polaron dispersion (the phonon being dispersionless) it is 
reasonable that as the polaron bandwidth decreases ex- 
ponentially with increasing coupling, so does the width 
of the polaron-|-one-phonon continuum. 

The higher energy features are also quite interesting. 
For the intermediate coupling A = 1.2, at some distance 
above the polaron-|-one-phonon continuum one can see 
the feature evolved from what was the third discrete state 
in the MA approximation. For a large 77 value this looks 
like a continuum, however scaling with 77 reveals a dis- 
crete state just below another continuum. In fact, the 
spectrum is broken up into discrete states and continua 
separated by gaps where no states are present. Of course, 
the detailed shape of the spectral weight above the third 
bound state is likely to change as one goes to MA(3) and 
higher orders, however it seems implausible that these 
gaps should all close and a single continuum should form 
above Eqs + as is the case at very weak couplings 
(see below). Instead, these results suggest that most 
weight is inside bound states which are reminiscent of 
the Lang-Firsov "comb" of discrete states separated by 
a frequency fi. Here the distance between discrete states 
is generally less than Q and there are narrow continua in 
between them, which evolve from lower-energy bound- 
states -|- one or more phonons. This is nicely illustrated 
by the results in Fig. [^b), which show 2 different fea- 
tures between the third and the fourth bound states. The 
lower-energy one is a continuum that starts roughly at f2 
above the second bound state, and has a finite width. 
Since the second bound state still has some finite band- 
width at this coupling (see below), it seems reasonable to 
interpret this feature as being the polaron in the second- 
bound state plus one phonon somewhere far from it. The 
higher energy feature is much narrower, but close inspec- 
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FIG. 7: (color online) (top) MA spectral weight A{0, lj) vs. 
ij in ID for t = 1, = 0.5, rj = 0.01 and A varying from 
to 2; (bottom) Same for MA'^'^-'. Curves corresponding to 
A = 0.3, 0.6, 1.2 and 1.8 are highlighted. 

tion reveals that it also is a continuum, with the lower 
edge starting at Eqs + 2ri, so its origin is obvious. The 
continua between higher bound states become wider, in 
agreement with expectations if one assumes that indeed 
they result from adding a distant phonon to a polaron 
in a higher bound state, which has a larger bandwidth. 
There is also overlap between different types of states, 
which also leads to increased bandwidth. 

This structure of the spectral weight explains why the 
MA (which predicts only bound states at low energies. 



for medium to strong coupling) still obeys sum rules with 
such good accuracy.^ Most of the weight is indeed in the 
bound states, not in the narrow continua that appear 
in between them. These results also show clearly how 
the convergence towards the Lang-Firsov limit g 3> i is 
achieved: the width and weight of the continua shrinks 
to zero, and one is left only with the equally spaced dis- 
crete states. This is illustrated in Fig. [3 where k = 
MA spectral weights are contrasted with MA'^' spectral 
weights for different A values. The largest difference is ob- 
served at small couplings, where MA predicts the (wrong) 
continuum pinned at —2dt + Q = —1.5 for these values, 
whereas MA^^^^ clearly shows a continuum starting at 
above the ground-state for as long as it is still visible. 
For A > 0.7 or so, a second bound state splits from this 
continuum and comes quite close to the GS peak before 
moving away so that it asymptotically goes to Eqs + ^■ 
There are clear similarities between the two plots, with 
most of the weight concentrated in the bound states that 
are roughly fl apart, but MA overestimates their weights 
in order to compensate for the missing narrow continua 
that appear in between these discrete states. As stated, 
the precise shape and weight of the higher continua is 
very likely to change as one goes to a higher level MA 
approximation, however we expect the general picture to 
remain essentially the same. 

We are unable to check quantitatively the accuracy of 
the higher energy spectral weight against detailed numer- 
ical predictions, beyond the sort of comparisons shown 
for the polaron dispersion in Fig. [S] The reason is that 
most of the numerical work is focused on computing low- 
energy properties (a detailed overview of such work is 
provided in Ref. |2, or, for example, in Ref. @). The 
much fewer high-energy results, ^ch as those based on 
a variational treatment in Ref. 10 and a novel QMC / 
exact diagonalization approach in Ref. [ill use a rather 
large ry and are already in reasonable agreement with 
MA, as discussed in Ref. 0. Cluster perturbation the- 
ory (CPT) results such as shown in Ref. [13, while also 
an approximation, are in good qualitative and quantita- 
tive agreement with ours. It would be very interesting 
to be able to compare our results against detailed high- 
accuracy, high-energy numerical predictions. 

Finally, we contrast the difference between MA and 
MA(2) spectral weights for different momenta and ener- 
gies. Typical results are shown in Figs. [5] and [21 for 
weak, intermediate and strong couplings A = 0.3,0.6, 1.2 
and 1.8, respectively. For the weak coupling, as already 
discussed, the most obvious difference is that the po- 
laron bandwidth is decreased to its correct value of 
m MA(2). The qp weights and all other features are 
very similar. Based on the MA'^^^ results, it is now clear 
that the strong resonance seen in the electron-|-phonon 
continuum occurs at an energy of 2U above the ground- 
state energy. One expects that here is where the second 
bound state will arise from. For A = 0.6, Figs. [5]Jc) and 
(d) show a bigger contrast. Here, MA already predicts 
a second bound state that has evolved in between the 
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FIG. 8: A{k, Lo) vs. fc and in ID for i = 1, f7 = 0.5, = 0.01 
and A — 0.3 in (a), (b) and A = 0.6 in (c), (d). Results for 
MA are shown in (a), (c), while MA*-^^ is shown in (b), (d). 
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FIG. 9: A{k, uj) vs. k and uj in ID iort = l,n = 0.5, = 0.01 
and A — 1.2 in (a), (b) and A — 1.8 in (c), (d). Results for 
MA are shown in (a), (c), while MA^^-* is shown in (b), (d). 



polaron band and the higher-energy continuum. In con- 
trast, MA*^^) shows that there is no second-bound state 
yet, however the electron-|-one-phonon continuum is spht 
off the higher energy continuum, which also starts to 
fractionaUze at higher energies (roughly multiples of Q). 
Within MA^^), a true second bound-state is observed for 
the higher couplings shown in Fig. [HI In contrast to 
MA, which shows several bound states which disperse 
as k increases, MA^^^ also clearly shows continua in be- 
tween these discrete states. These account for some of the 
spectral weight that was in the MA peaks. These results 
again suggest a very fractionalized spectrum at interme- 
diate and strong coupHngs. Instead of the polaron band, 
a second-bound state and a rather featureless continuum 
at higher energies, we instead find many sets of discrete 
states interspersed with continua. As A ^ oo, the rela- 
tive weight in these continua decreases and the spectral 
weights evolve continuously towards the Lang-Firsov set 
of discrete states with energies —g^/il + nft. 

Qualitatively similar results are observed in higher di- 
mensions. Rather ironically, the most time-consuming 
part in the MA*^^^ calculation for higher dimensions is 
finding the spatial dependence of the non-interacting 
Green's functions Go{i,n,) which are needed to generate 
Eqs. P5)) . and not the solving of the system of equations. 
The reason is that for nearest-neighbor hopping in higher 
dimensions, the evaluation of these propagators must be 



done by numerically. Of course, one could choose simpler 
forms of the dispersion for which analytical results are 
possible. However, as we show now, in higher dimen- 
sions the improvements in going to MA*^^) and MA^^^ are 
quantitatively smaller, because the relative weight in the 
continua is reduced compared to the ID case. This is in 
agreement with our general observation that MA itself 
becomes more accurate with increased dimensionality!^ 

In Figs. [T0lfT3l we show the 2D spectral weight v4(k = 
0,uj) vs. UJ, for effective couphngs A = 0.3,0.6,0.95 and 
1.2, both on linear and logarithmic scales. Qualitatively, 
everything is similar to the behavior seen in the ID case. 
Quantitatively, we find that the continuum that appears 
at fl above the ground-state is broader, but of lower 
height. In fact, its height is so small that it is invisible 
on curves like those in Fig. [71 which is the reason why 
we do not show such curves here. The overall weight in 
this continuum also decrease much faster with increasing 
A. In ID, for A = 1.2 the first contimmm is still clearly 
visible, see Fig. O For A = 1.8 it becomes harder to see 
on the linear scale, but it is clearly seen on the logarith- 
mic scale. By contrast, in 2D, for A ~ 1.2 the continuum 
is no longer visible on the linear scale, and even in the 
logarithmic scale it only barely starts to be visible for 
77 = lO""*^ (small shoulder marked by arrow). The spec- 
tral weight for a Lorentzian contribution Z/(lu + irf) is 
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FIG. 10: (color online) (a) A{Q,u)) vs. uj in 2D for t = 
l,n = 0.5, A = 0.3,77 = 0.01, in MA, MA^^' and MA'^' 
(curves shifted for clarity); (b) lnyl(0,Ci;) vs. lo for MA(2) 
and r\ = 10~^, 10~^, 10~^. Other parameters are as in (a). 



FIG. 12: (color online) (a) A(0,a;) vs. u in 2D for f = 1, H = 
0.5, A = 0.95,7? = 0.01, in MA, MA'^' and MA^^' ( curves 
shifted for clarity); (b) lnyl(0, u) vs. u for MA'^> and r? = 
10~^, 10"'', 10"*. Other parameters are as in (a). 



Zr\l\Tx(LL)^ + so the maximum height of the peak, at 
resonance, is Z /■nrj. This is visible only if it is larger than 
the background, which of course depends on how close is 
the next spectral feature. However, for a very small Z , i] 
has to be really small before the peak is seen. 

As in the ID case, we also observe the fractionaliza- 
tion of the spectrum for moderate and large effective 
couplings, with a succession of discrete peaks and con- 
tinua at higher energies. We expect similar behavior to 
be observed in 3D as well. Clearly, the changes in going 
from MA to MA^^' and MA^^^ are quantitatively much 
smaller in 2D than in ID, because the continua have so 
little weight. We expect the trend to continue in going 
to 3D, meaning that in 3D, the difference between MA*^^^ 
and MA should be quantitatively even less. Indeed, all 
the comparisons of MA results with available 3D numer- 
ics, shown in Ref. are already in excellent agreement, 
even for intermediary couplings A 1. As a result, we 
find it unnecessary to present 3D results, although they 
can be obtained very straightforwardly. 



V. CONCLUSIONS 

In summary, we presented a way to systematically im- 
prove the MA approximation, by systematically improv- 
ing the accuracy of self-energy diagrams, but in such a 
way that they can still all be efficiently summed. 

This allows us to rather easily fix various known fail- 
ings of the MA approximation, such as the absence of 
the polaron-|-one-phonon continuum at the correct en- 
ergy, and its momentum-independent self-energy. It also 
allows us to understand in more detail the effects of 
the Holstein-type electron-phonon coupling on the po- 
laron spectrum, both at low and high energies. While 
agreement with exact numerical results is improved, un- 
fortunately there are not many such results for higher 
energy states, so detailed comparisons are not possible 
there. However, the hierarchy of MA^") approximations 
is clearly providing a simple way towards obtaining quan- 
titatively more and more accurate results for the Green's 
function of the Holstein polaron, in any dimension and 
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FIG. 11: (color online) (a) A(0,l<;) vs. in 2D for t = 
1,^ = 0.5, A = 0.6,77 = 0.01, in MA, MA^^' and MA*^' 
(curves shifted for clarity); (b) lnyl(0,a;) vs. lo for MA(2) 
and 77 = 10~^, 10~^, 10~*. Other parameters are as in (a). 
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FIG. 13: (color online) (a) A(0,a;) vs. in 2D for t = 
1,J7 = 0.5, A = 1.2,77 = 0.01, in MA, MA^^' and MA*^^ 
(curves shifted for clarity); (b) lnj4(0,t<;) vs. 10 for MA(2) 
and r\ — 10~^, 10~^, 10"*. Other parameters axe as in (a). 
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for any free electron dispersion. 

The next direction of obvious interest is to study to 
what extent this work can be extended to other models, 
for example models where the electron-phonon coupling 
depends on the phonon momentum {e.g. the breathing- 
modei^ or Frolicbi^ Hamiltonians). For such models 
there are very few reliable high-energy results available. 
A simple approximations such as MA could easily inves- 
tigate the whole parameter space and identify interesting 
regimes. Such work is currently in progress. 
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APPENDIX A: CONTINUED FRACTIONS 

Consider the set of recurrence relations /_i = 1 and 
for n > 0, 



/n — CKn/n-l + Ai/n+l, 

where /„ — s- as n — s- oo (this is justified for all cases of 
interest to us since /„ is related to a generalized Green's 
function with n phonons in the initial state and none in 
the final state. If n is large enough, the amplitude of 
probability to evolve between such states must vanish). 

Truncating the system at a large enough N so that 
/n+i ~ 0, we solve backwards to find /at = aNfN-i, 
/n-i ^ aN-ifN-2/{l - Pn-i^n), In general. 



^nfn- 



1 - 



1 - 



1 - 



If we allow N 

An = a„/[l - 



— > oo, then the solution becomes exact and 
[3nAn+i\ are infinite continued fractions. 



so that Sjv^^(i) (k, w) = J-i. 

First, we average each of Eqs. pS]) over all the cor- 
responding momenta, to find that for n > 2, Tn — 
go{uj — nfl)[ng^J^n-i + ^n+i]- This has the continued- 
fraction solution (see Appendix E]) : 



(B3) 



APPENDIX B: MA^^^ DETAILS 



g^A2iiu)Ti 



where A2{uj) is a continued fraction defined in Eq. ([l3 
and J^i is still unknown. 

We now average each of Eqs. over all the momenta 
except qi, to find that for all n > 2, /n(qi) — go{uj — 
nn)[g^J^n-i + {n-l)g'^fn-i{q.i)+fn+i{<li)]- This can also 
be solved by rewriting it in terms of (5/„(qi) = /n(qi) — 
JF„. Using the recurrence relation for J^n, we find that: 

Sfn{q.i) = .go(^ - nl7)[(n - l).g^J/„_i(qi) -|- 5/„+i(qi)]. 

This also has solutions in terms of continued fractions, 
in particular (5/2 (qi) — g^Ai{uj — r2)(5/i(qi). Using 
Eq. dm) we then find /2(qi) = g^Ai{uj - rj)/i(qi) + 
g^ [A2{iv) — Ai{u; — JFi. This expression is now in- 
serted in Eq. p4)) to give an equation with only /i"^^(qi) 

and Ti — J2qi fi^\^i) unknowns, which can be 
re- arranged as: 



(Go(k - qi,c. - n))-' - g^A.iuo - VL)] ^(qi) 
= 5'+ff' [A2{uo)~A^{uo^9)]Ti. 



But [Go(k,a;)] ^ — a{ijj) = oj — ek + irj — a{u!) = [Go(k, 
a{uj))]~^ and therefore we find: 

/^'^(qi) = Go(k - qi, - r! - g^A^ioj - il)) 



X [g^+g^ [A2iLj)~AiiLj^n)]J^i] . 

After momentum averaging over qi on both sides, we 
find JFi = S [see Eq. ([T5)) ]. Based on the knowledge of 
f[^\cii) and the various partial and total averages, one 
can now compute /2^''(qi, q2), etc. 



First, we use Eqs. ([TS]) to find how J2q2 /2^^(qi' ^2) 

depends on f[^\c[i). This is then used in Eq. ^4]) to 

solve for /i^''(qi), and S follows from Eq. p^ . 
We introduce the partial momentum averages: 

/n(qi)-^ E /rl'Hqi,---,qn), (Bi) 

q2,..,q,i 

(we use the convention that /i(qi) — /j^^(qi)) and the 
full momentum averages: 

qi,--,qn qi 

(B2) 



APPENDIX C: MA'^' DETAILS 

As for MA(^\ the strategy is to solve Eqs. to 

find how ^ /i^^qi' q2: qs) depends on /^^^(qi, q2). 
Using this in Eq. (|20p allows us to solve it in conjunction 

with Eq. HH) to find f['^\qi) and therefore S. 
We introduce various partial momentum averages: 

/n(qi,q2) = /,l^Hqi,q2,-.-,qn), 
q3,--,qn 

/n(qi) = fnH^lU ■ ■ • ,qn)- 

q2,--,q,i 
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where we take /i(qi) = /i(qi) = /} (qi),/2(qi, q2) = 

(2) 

f2 Qs)- We also define the full momentum averages: 

•^" = 4 E /l''(qi,---,q«) = ^E/"(qi) (CI) 



The functions /„,/n and ^„ should also carry the up- 
per label (2) since their values are different than those 
obtained within MA^^^ however we do not write it ex- 
plicitly to simplify the notation somewhat. 

The solutions for J-„ and /„ proceed just as in the 
case, with one difference. Averaging all equations 
()2ip over all corresponding momenta, we find that now 
only for n > 3, Tn — gQ{Lu — nQ,)[ng'^Tn-i + J^n+i], there- 
fore this recurrence relation now ends with: 



(C2) 



Similarly, after using the same approach as for MA^^^ 
one jinds that 5/3(q_i) = /sCqi) - -^3 = g^^2('^ - 
r!)<5/2(qi) so that /3(qi) - g^A2{LJ - r!)/2(qi) + 
g'^ [A3{u;) — A2{u! ~ n)] ^2- Of course, the functions 
J-2,/2(qi) are still unknown. Finally, we can proceed 
to calculate /3(qi,q2) which is needed in Eq. ([20|l . 

We momentum average Eqs. ([2T|) over all momenta 
except qi , q2 , to find 

/«(qi,q2) = goiuJ - nn) [g^ (/„-i(qi) + /„-i(q2)) 
2).g2/„_i(qi,q2) + /„+i(qi,q2) • 

Let 6fn (qi , q2 ) = (qi , q2 ) - In (qi )- In (q2 )+^n- For 
these, the recurrence relations are 5/n(qi,q2) = go{^ — 
nn) (n - 2)g2(5/„_i(qi,q2) + (5/„+i(qi,q2) with the 
solution <5/3(qi^q2) = g^Ai{uj - 2n)6l2{qi,c[2), from 
which we find /3_(qi,q2) = ff^Ai/2^^(qi, q2) + g'^{A2 - 
Ai) [Shini) + (5/2 (q2)] + <?' [^3 - ^i] ^2, i.e. only 

in terms of various partial averages of (qi,q2). 
Throughout we use the shorthand notation 

Ai = Aiiuj - 2n),A2 = A2iuj - n),A3 = Asiuj). 
We now substitute /3(qi,q2) in Eq. ((20|) . which gives: 

/2(qi,q2) - g'Go(k- qi - q2,cD) [/f ^(qi) + /f ^(q2) 
+{A2 - Ai) [(5/2(qi) + 5/2(q2)] + [As - A,] T2\ . (C3) 

We used the fact that [Go(k— qi — q2, 2f2)]~^ — g^Ai — 
[Go(k — qi — q2, "i)]""^, with d; = w — 2f2 — g^ A\. 



Since only /2(qi) is needed in Eq. p^ . we can obtain 
it by momentum averaging Eq. (jC3[) over q2 . This leads 
directly to the closed system of coupled equations: 

/f ^(qi) = Go(k - qi,c. - f^) [g2 + ^/^(qi) + , 



^/2(qi) = .9'.9o('^) /r'(qi) + (^2 - ^i)'5/2(qi) ~ 



' N 



^Go(k-qi-q2,(:>) \f['\q2) + {A2-Ai)6l2{q2) 



To solve it, we first rewrite these as a single equation in 
terms of the unknowns Xq = /{ (q) -I- {A2 — Ai)Sf2{q)- 
This is obtained by adding the second equation to [A2 — 
Ai + Go(k — qi,uj — O)] times the first equation. We 
introduce the short-hand notation: 

fly = aij{uj) = 1 - g^go(i^)[Ai - Aj] 

and use the identities [Go(k — qijW — f^)]^^ ^ 
g2go(ci)[l + (A2-Ai)(Go(k-qi,c^-f^))-i] = [1 - 

5r2go((I,)(A2 - ^i)] [Go(k-qi,(I>)] , where 



a2i 



and Go(k-qi,(5) 1 -I- (A2 - Ai) [Go(k - qi, w - fl)]" 

A2 — Ai + Go(k — qi,(D)/a2i. Both identities are based 
on the definition [Go(k, uj)]~^ — to — + irj. 
With these, the final equation for Xq becomes: 



a2iXq = Go(k- qi,^) 



-031(^2 -Ai)J^2 



A2-A1 



2 a2i — 031 

g H T2 

021 

Go(k- qi,(I;)' 



a2i 



N 



^Go(k- q - q2,<I')a 



12 ■ 



q2 



We now Fourier transform x{i) = ;^ ^''"^'^q- 
From the definition of Xq it follows immediately that 
x{0) = = E. Since J^2 ^1 = x{0) [see Eq. (IC4l) ]. 
the resulting system of inhomogeneous equations is: 



where Go{i,uj) — e*'^'^'Go(k, a;), and the matrix 

elements are: 



Momentum averaging Eq. (|C3P over both momenta 
leads immediately to JF2 = 5^3o(w) [2J^i + (A3 — Ai)J^2], 
and therefore: 



Moo = 1 - g^ffo('^)5o('^ 
M»o = -g'5o('^)e''"''*GoH,(i) 



' 2 


^) 




a2i / 


(- 


1 ^ 




021, 



(C5) 
(C6) 



^2 



2g^go{Ci^)Ti 



l-.g2.go((:;)[yl3-yli] 



(C4) 



It follows that all higher order total averages are pro- 
portional to J- 1 — Ti. 



for i ^ 0, and for both i,jy^O 

My =a2i<5.,, -ffV^-^'Goad-) 

{A2 - Ai)6i^_j + 



2„ikRi 

Go{~i~ j,oj) 



0,21 



(C7) 
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We truncate this system by keeping only sites R; within 
a certain distance from the origin (see discussion in main 
text) and solve it numerically. The self-energy is 

5^MA(2)(k,w) = X'(0). 
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